#include "header.h"

#define NONB    10000
#define MINRAT	1.0e-4

double MaxSVD, MinSVD;
/* ----------------------------------------------------- */
void svd_test(struct node* p){
/* ----------------------------------------------------- */
int dsvdcmp(double **a, int m, int n, double *w, double **v);
void FindGeoFull(int, struct node*, double**, double*, double**);
void challoc(void*);

int i,j,neigh, status;
double dx,dy,dz,weight,ds2;
double **AF,w[3+1],**v;
FILE* fpt;

printf("Performing SVD test, see SVD.OUT ...\n");

AF  = (double**)calloc(NONB,sizeof(double*));
MinSVD = 1.0e20;
MaxSVD = 0.0;
for(j=0;j<NONB;j++){
	AF[j]  = (double*)calloc(3+1,sizeof(double));
	}
v = (double**)calloc(3+1,sizeof(double*));
for(j=0;j<3+1;j++) v[j] = (double*)calloc(3+1,sizeof(double));

fpt = fopen("SVD.OUT", "w");
for(i=0;i<no_grid_points;i++){
	if(p[i].total_conn > NONB){
		printf("svd_test: Too many neighbours. Increase size of NONB.\n");
		exit(0);
		}
	for(j=0;j<p[i].total_conn;j++){
		neigh = p[i].Conn[j];
		dx    = p[neigh].x - p[i].x;
		dy    = p[neigh].y - p[i].y;
		dz    = p[neigh].z - p[i].z;
		ds2   = dx*dx + dy*dy + dz*dz;
		/* Note that difference in the power */
		weight= sqrt( WEIGHT(ds2, p[i].ds2) );
        //weight = 1.0;

		dx    = weight*dx;
		dy    = weight*dy;
		dz    = weight*dz;

		if(p[i].point_type == SOLID_WALL || //p[i].point_type == OUTER_BOUND ||
		   p[i].point_type == FIRST_LAYER)
			rotate_x(&dx,&dy,&dz,&p[i]);

		AF[j+1][1] = dx; AF[j+1][2] = dy; AF[j+1][3] = dz;

		}


	/* Coefficients for Full stencil */
	status = dsvdcmp(AF, p[i].total_conn, 3, w, v);
    if(status == FAILURE){
         printf("Full stencil\n");
         printf("Node = %d, conn = %d, type = %d\n", i, p[i].total_conn, 
                p[i].point_type);
         fflush(NULL);
         exit(0);
         }
	FindGeoFull(i, p, AF, w, v);

    fprintf(fpt, "%d %e %e %e\n", i, w[1], w[2], w[3]);
	}

fclose(fpt);

/* Free the allocated memory */
for(j=0;j<NONB;j++){
	free( AF[j]  );
	}
free(AF);
for(j=0;j<4;j++) free(v[j]);
free(v);

printf("Minimum svd = %e\n",MinSVD);
printf("Maximum svd = %e\n",MaxSVD);
printf("MinSVD/MaxSVD = %e\n",MinSVD/MaxSVD);

/* xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx */
}
/* -------------------------------------------------------------------------- */
void FindGeoFull(int node, struct node* p, double** U, 
             double* W, double** V){
/* -------------------------------------------------------------------------- */
double wmin, wmax;
if( W[1] == 0 || W[2] == 0.0 || W[3] == 0.0){ /* Required derivative cannot be determined */
	printf("svd_test: Fatal error, full stencil insufficient\n");
	printf("           for determining derivatives at node %d.\n",node);
	printf("           x= %e\n", p[node].x);
	printf("           y= %e\n", p[node].y);
	printf("           z= %e\n", p[node].z);
	printf("           Point type= %d\n", p[node].point_type);
	printf("           Conn size = %d\n", p[node].total_conn);
	exit(0);
	}

if( W[1] < 0 || W[2] < 0.0 || W[3] < 0.0){ /* Required derivative cannot be determined */
	printf("FindGeoFull: Fatal error, svd is negative at node %d\n",node);
	printf("             Point type = %d\n",p[node].point_type);
	exit(0);
	}

wmin = min(W[1], W[2]);
wmin = min(wmin, W[3]);

wmax = max(W[1], W[2]);
wmax = max(wmax, W[3]);

MinSVD = min(MinSVD, wmin);
MaxSVD = max(MaxSVD, wmax);

if(wmin/wmax < 1.0e-3){
	printf("FindGeoFull: FATAL error at node %d\n", node);
	printf("             Full stencil is ill-conditioned\n");
	printf("             %e %e %e\n", W[1], W[2], W[3]);
	//exit(0);
	}

}
/* ------------------------------------------------------------ */
void challoc(void* var){
/* ------------------------------------------------------------ */
	if(var==NULL){
		printf("Memory allocation failure.\n");
		exit(0);
	}
}
/* ------------------------------------------------------------ */
